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Abstract 



Low-order equal-time statistics of a barotropic flow on a rotating sphere are 
investigated. The flow is driven by linear relaxation toward an unstable zonal jet. 
For relatively short relaxation times, the flow is dominated by critical-layer waves. 
For sufficiently long relaxation times, the flow is turbulent. Statistics obtained from 
a second-order cumulant expansion are compared to those accumulated in direct 
numerical simulations, revealing the strengths and limitations of the expansion for 
different relaxation times. 
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1. Introduction 



Many geophysical flows are subject to the effects of planetary rotation and to forcing and 
dissipation on large scales. For example, the kinetic energy of atmospheric macroturbulence is 
generated by baroclinic instability and is then partially transferred to mean flows, whose energy 
dissipation can often be represented by large-scale dissipation. Statistically steady states of such 
flows can exhibit regions of strong mixing that are clearly separated from regions of weak or no 
mixing, implying that the mixing is non-ergodic in the sense that flow states are not phase- space 
filling on phase space surfaces of constant inviscid invariants such as energy and enstrophy 
( Shepherd] | 1987| ). As a consequence, concepts from equilibrium statistical mechanics, which 
rely on ergodicity assumptions and can account for the statistics of two-dimensional flows in 



the absence of large-scale forcing and dissipation (e.g., Miller 1990 Robert and Sommeria 



1991[|Turk ingto n~et al.||2001||Majda and Wang 2006), generally cannot be used in developing 



statistical closures for such flows. 

In this paper, we investigate the inhomogeneous statistics of what may be the simplest flow 
subject to rotation, large-scale forcing, and dissipation that exhibits mixing and no-mixing re- 
gions in statistically steady states: barotropic flow on a rotating sphere driven by linear relax- 
ation toward an unstable zonal jet. Depending on a single control parameter, the relaxation time, 
this prototype flow exhibits behavior in the mixing region near the jet center that ranges from 
critical-layer waves at short relaxation times to turbulence at sufficiently long relaxation times. 
This permits systematic tests of non-equilibrium statistical closures in flow regimes ranging 
from weakly to strongly nonlinear. 

We study a non-equilibrium statistical closure based on a second-order cumulant expansion 
(CE) of the equal-time statistics of the flow. The CE is closed by constraining the third and 
higher cumulants to vanish, and the resulting second-order cumulant equations are solved nu- 
merically. The CE is weakly nonlinear in that nonlinear eddy-eddy interactions are assumed to 
vanish. We show that for short relaxation times, the CE accurately reproduces equal-time statis- 
tics obtained by direct numerical simulation (DNS). For long relaxation times, the CE does not 
quantitatively reproduce the DNS statistics but still provides information, for example, on the 
location of the boundary between the mixing and the no-mixing region — information that local 
closures (e.g., diffusive closures) do not easily provide. 

Section [2] introduces the equations of motion for the flow and discusses their symmetries 
and conservation laws. Section [3] describes the DNS, including the accumulation of low-order 
equal-time statistics during the course of the simulation. The CE and its associated closure 
approximation are outlined in section |4j Section [5] compares DNS and CE. Implications of the 
results are discussed in section |6] 



2. Barotropic jet on a rotating sphere 

a. Equations of motion 

We study forced-dissipative barotropic flow on a sphere of radius a rotating with angular 
velocity Vt. Though not crucial for this paper, we prefer to work on the sphere and not in the 
/5-plane approximation, as the sphere can support interesting phenomena not found on the plane 
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(e.g., |Cho and Polvani|1996| ). The absolute vorticity q is given by 



q = C + f 

= V 2 ^ + / (1) 
where £ is the relative vorticity, ip is the stream function, V 2 is the Laplacian on the sphere, and 

/(0) = 2fisin0 (2) 

is the Coriolis parameter, which varies with latitude 0. The time evolution of the absolute 
vorticity is governed by the equation of motion (EOM) 

— + J[i/;,q] = ^— — , (3) 

where 

j\ I ]- 1 fdipdq dipdg\ 

[W) ql ~ a 2 cos(0) V d\ <90 d<f> d\ J W 

is the Jacobian on the sphere and A is longitude. Forcing and dissipation are represented by 
the term on the right-hand-side of Eq. ([3]), which linearly relaxes the absolute vorticity q to the 
absolute vorticity q- ]Ct of a zonal jet on a relaxation time r. 

The zonal jet is symmetric about the equator and is characterized by constant relative vor- 
ticities ±T on the flanks far away from the apex and by a rounding width A0 of the apex, 

g jet (0) = /(0)-rtanh(^ . (5) 

The meridional profile ([5]) in our simulations is shown in Fig. [5] below. In the limiting zero- 
width case A0 — > of a point jet, 

Cjct(0) = g jot (0) - /(0) = -r sg n(0) , (6) 

and the jet velocity has zonal and meridional components 

Mjet (0) = ratan(|0|/2-vr/4), 

M0) = 0. (7) 

For T > 0, the zonal velocity attains its most negative value — Ta at the equator. 

For T > 0, the gradient of the absolute vorticity ([5]) changes sign at the equator, so the jet 



satisfies the Rayleigh-Kuo necessary condition for inviscid barotropic instability. Lindzen et al. 
(1983) showed that the linear stability problem for the barotropic point jet on a /3-plane is ho- 
momorphic to the Charney problem for baroclinic instability. The analog of the horizontal sign 
change of the absolute vorticity gradient in the barotropic problem is the vertical sign change 
of the generalized potential vorticity gradient in the baroclinic problem (with the generalized 
potential vorticity gradient including a singular surface contribution from the surface potential 
temperature gradient). What is the meridional coordinate in the barotropic problem corresponds 
to the vertical coordinate rescaled by the Prandtl ratio in the baroclinic problem. 
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The analogy to the baroclinic Charney problem motivated extensive study of the barotropic 
point-jet instability and its nonlinear equilibration, with the forcing and dissipation on the right- 
hand side of Eq. ([3]) as a barotropic analog of radiative forcing by Newtonian relaxation and 
Rayleigh drag (e.g., |Schoeberl and Lindzen||1984| [Nielsen and Schoeber l|1984[ [Schoeberl and 



Nielsen 1986; Shepherd 1988). Building on this body of work, here we focus on the statistically 



steady states of the flow and study their dependence on the relaxation time r. This allows us 
to test non-equilibrium closures systematically in two-dimensional barotropic flows that exhibit 
similar phenomena as analogous three-dimensional baroclinic flows, with the caveat, of course, 
that additional degrees of freedom in three-dimensional baroclinic flows, such as adjustment 
of the static stability, make the equilibration of baroclinic instabilities different from that of 
barotropic instabilities. 

b. Symmetries and conservation laws 

Because the jet to which the flow relaxes is symmetric about the equator, steady-state statis- 
tics of the flow are hemispherically symmetric. Deviations from hemispheric symmetry can 
be used to gauge the degree of convergence towards statistically steady states. They will also 
highlight a qualitative problem with the statistics calculated by the CE (see section [5] below). 

The EOM, Eq. ([3]), is invariant under a rotation of the azimuth, A — > A + a, and under an 
inversion of the coordinates, 



A -> -A 
q -»• q 

<?jet -> ?jet • (8) 

Furthermore, the vorticities change sign under a north-south reflection about the equator, 

A -> A 

q -»• -q 

<?jet -> -<?jet • (9) 

These symmetries are reflected in the statistics discussed below. 

As a consequence of the constancy of the relaxation time r, statistically steady states satisfy 
two constraints that can be obtained by integrating the EOM over the domain. Kelvin's circula- 
tion and Kelvin's impulse of long-time averages (■) in a statistically steady state are both equal 
to those of the jet to which the flow relaxes, 



(q(r,t)) dr = J q jet dr, (10) 
J (q(r,t)) sin0 dr = J g jct sin <\> dr, (11) 



where r = ((f), X) is a position vector. Conservation of circulation (10) is trivially satisfied 
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because vorticity integrals vanish at each moment in time, 



J q(r,t) dr = J C(r,t) dr = J q iet (r,t) dr = 0. (12) 

However, conservation of impulse (jTT]), which is equivalent to conservation of the angular mo- 
mentum about the rotation axis, is not trivial and must be respected by statistical closures. 

3. Direct numerical simulation 

a. Parameters and implementation 

All vorticities and their statistics can be expressed in units of VL, but to give a sense of 



scale, we set the rotation period to 2-k/VL = 1 day. We use Arakawa's (1966) energy- and 
enstrophy-conserving discretization scheme for the Jacobian on a M x N grid. For all results 
reported below, there are M = 400 zonal points and N = 200 meridional points. The lattice 
points are evenly spaced in latitude and longitude, apart from two polar caps that eliminate the 
coordinate singularities at the poles. Each cap subtends 0.15 radians (8.6°) in angular radius 
and, following Arakawa, consists of a union of triangles radiating from the pole that match the 
grid along their base; scalar fields are constrained to be constant in each cap. At the initial time 
t = 0, we set q = q- ]et plus a small perturbation that breaks the azimuthal symmetry and triggers 
the instability. The time integration is then carried out with a standard second-order leapfrog 
algorithm using a time step of At = 15 s. The accuracy of the numerical calculation was 
checked, in the absence of the jet, against exact analytic solutions that are available for special 



initial conditions ( |Gates and Riegel|1962[ ). The jet parameters are fixed to be T = 0.6Q and 



A<fi = 0.05 radians (2.9°). Though unphysically fast for Earth, the jet illustrates the strengths 
and shortcomings of the CE. Code implementing the numerical calculation is written in the 
Objective-C programming language, as its object orientation and dynamic typing are well suited 
for carrying out a comparison between DNS and the CE. 

Figure [T] shows the zonal velocity of the fixed jet «j e t(0) an d the mean zonal velocity of the 
flow (u((fi)) for two different relaxation times. Rounding of the jet due to mixing is evident. The 
absolute vorticity during the evolution of the instability and in the statistically steady state even- 
tually reached in a typical DNS are shown in Fig. [2j Figure|3]displays snapshots of the absolute 
vorticity in the steady-state regime for six different choices of r. In the limit of vanishingly 
short relaxation time r — > and strong coupling to the underlying jet, the fixed jet dominates, 
and q = gj et with no fluctuations in the flow. For r > 0, instabilities develop, and irreversible 
mixing begins to occur in critical-layer waves, which form Kelvin cats' eyes that are advected 
zonally with the local mean zonal flow (e.g., |Stewartson||1981||Maslowe||1986| ). At sufficiently 



large relaxation times (r > 12 days), the jet becomes turbulent, and as r increases further, tur- 
bulence increasingly homogenizes the absolute vorticity in a mixing region in the center of the 
jet. The dynamics are strongly out of equilibrium and nonlinear for intermediate values of r, 
yet continue to be statistically steady at long times. In the limit of long relaxation time r — > oo 
and weak coupling to the underlying jet, and upon addition of some small viscosity to the EOM, 
the system reaches an equilibrium configuration at long times (Salmon 1998; Turkington et al.| 



2001 1 Weichman 2006, |Majda and Wang 2006), and again the fluctuations vanish. Here we 



restrict attention to the geophysically most relevant case of short and intermediate relaxation 
times. 

Part of what makes this flow an interesting prototype problem to test statistical closures is 
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that, except in the extreme limits of vanishing or infinite relaxation time, irreversible mixing is 
confined to the center of the jet and does not cover the domain. An estimate of the extent of 
the mixing region can be obtained by considering the state that would result by mixing absolute 
vorticity in the center of the jet such that it is, in the mean, homogenized there and continuous 
with the unmodified absolute vorticity of the underlying jet at the boundaries of the mixing 
region. Because of the symmetry of the jet, this state would have mean absolute vorticity 

f for 101 < <b c , 

9 H f U\Z7 (13) 

[q }ct for |0| > <p e , 

and the boundaries of the mixing region would be at the latitudes at which q^ ct = 0, which are, 



with our parameter values, C r/(20) ~ 17° (cf. Schoeberl and Lindzen 1984, Shepherd 
1988 )Q The meridional gradient of the resulting mean absolute vorticity does not change sign, 



so the corresponding flow would be stable according to the Rayleigh-Kuo criterion. It represents 
a zonal jet that is parabolic near the equator. However, while the mean absolute vorticity satisfies 



the circulation constraint (10) not only in the domain integral but integrated over the mixing 



region between ±0 C , it does not satisfy the impulse constraint (11). To satisfy the impulse 
constraint, the mixing region in a statistically steady state extends beyond the latitudes C , as 
can be seen in Fig.|3]and will be discussed further below. Statistical closures must account for 
the structure of the transition between the mixing and no-mixing regions in this flow. 

b. Low-order equal-time statistics 

The first cumulant (or first moment) c\ of the relative vorticity depends only on latitude 0, 
reflecting the azimuthal symmetry of the EOM, 

Cl (r) = (C(r)) = Cl (0). (14) 

It is also convenient to define the first moment of the absolute vorticity, qi(4>) = (q(r)) = 
Ci(0) + f{4>). The calculation of the time averages (•) commences once the jet has reached a 
statistically steady state. As the adjustment of the mean flow is controlled by the relaxation time 
r, reaching a statistically steady state takes longer for larger r. Statistics are then accumulated 
every 100 minutes for a minimum of 100 days of model time, until adequate convergence is 
obtained. We have verified that the long-time averages thus obtained are typically independent 
of the particular choice of initial condition; see, for instance, Fig. |4j (The one exception is 
the case of r = 3.125 days, in which, depending on initial conditions, critical-layer waves of 
wavenumber three or four are present in the statistically steady states. As the wavenumber- 
three mode has slightly lower kinetic energy, and is reached from strongly perturbed initial 
conditions, we focus on it here.) As expected, azimuthal symmetry is recovered in such long 
time- averages, as can be seen, for instance, in the final panel of Fig. [2j In addition, the first 
cumulant changes sign under reflections about the equator, 

Cl (-0) = - Cl (0) , (15) 

a consequence of the reflection symmetry Q. 



'in the analogy to the Charney problem, the meridional scale oT/(2f2) = T/j3 with j3 = 2fl/a is the barotropic 



analog of the vertical Held (1978) scale over which quasigeostrophic potential vorticity fluxes associated with 



unstable baroclinic waves extend. 
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The second cumulant of the relative vorticity, given in terms of its first and second moments 

by 



c 2 (r, r') = (C(r) C(r'))c = (CW C(r')> " (CW> (C(r')> , (16) 
depends on the latitude of both points r and r', but only on the difference in the longitudes: 

c 2 (r,r') = 02(0,0', A -A'). (17) 



It is essential to take advantage of the azimuthal symmetry of the second cumulant, Eq. (17), 
to reduce the amount of memory required to store the second cumulants by a factor of M, from 
M 2 N 2 to MN 2 scalars. In the DNS, the reduction is realized by averaging the second cumulant 
over A' for each value of A A = A — A'. The averaging also improves the accuracy of the statistic. 

By definition, the second cumulant is symmetric under an interchange of coordinates, c 2 (r, r') 
C2(r', r). It also possesses the discrete symmetry 

c 2 (-0, -0', AA) = c 2 (0, 0', AA) , (18) 

a consequence of the reflection operation ([9]). 

4. Second-order cumulant expansion 

The jets considered here are inhomogeneous and possess nontrivial mean flows. As a con- 
sequence, the leading-order nonlinearity is the interaction between the mean flow and eddies 
(fluctuations about the mean flow), and already the mean flow, a first-order statistic, is of inter- 
est in closure theories (e.g. |Schoeberl and Lindzen|1984[ ). In contrast, the leading-order nonlin- 
earity in homogeneous flows is the interaction of eddies with each other, and only higher-order 
statistics are of interest in closure theories. Many higher-order closure approximations impose a 
requirement of homogeneity (e.g., |Holloway and Hendershott||1977| |Legras||1980[ |Huang et"al 



2001 ) and are not directly applicable to systems with inhomogeneous flows. 



The second-order closure we consider is based on an expansion of vorticity statistics in 
equal-time cumulants. The cumulant expansion can be formulated either using the Hopf func- 
tional approach ( Frisch|1995 Ma and Marston|2005 ) or by a standard Reynolds decomposition 



of each scalar field into a mean component plus fluctuations or eddies (denoted with a prime): 

q(r) = ( C (r)) + /(0) + C'(r) 

V>(r) = (V>(r)) + V»'(r) ■ (19) 

The EOMs for the first and second cumulants may be written most conveniently by introducing 
the following auxiliary statistical quantities: 

Pi(r) = (V>(r)) 
p 2 (r,r') = (V(r)C(V)>c 

= WW). (20) 

These quantities contain no new information as c\ = V 2 p\ and c 2 = V 2 j?2 ? where it is un- 
derstood that unprimed differential operators such as V 2 and J[ , ] act only on the unprimed 
coordinates r. Substituting Eq. ( 19 ) into Eq. ^ and applying the averaging operation (•) yields 
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the EOM for the first moment or cumulant: 



dCl(r) ~ i[ci(r) + /W lPl (r)] + /^(r-r'),ft(r 1 r'Pr'+ (ictW ° i{t) 

(21) 



Here partial integration over r' has been used to group and (' that appear as separate ar- 
guments of the Jacobian operator into the statistical quantity p 2 . Similarly, multiplication of 
Eq. by C'(r') followed by averaging yields the EOM for the second cumulant, which upon 
discarding the term cubic in the fluctuations may be written as 

dC2 ^ r>) = J[c 1 (r) + /(0),p 2 (r,r')] + J[c 2 (r,r / ), Pl (r)] - + (r ~ r'), 

(22) 

where (r <-» r') is shorthand notation for terms that maintain the symmetry 02(1", r') = c 2 (r', r). 
Closure has been achieved at second order in the expansion by constraining the third and higher 
cumulants to be zero 

c 3 = (C(r) C(r') C(r"))c = 0, etc. (23) 

Otherwise an additional term would appear in Eq. ( [22] ) that couples the second and third cumu- 
lants. The closure approximation c 3 = amounts to neglecting eddy-eddy interactions while 
retaining eddy-mean flow interactions (e.g., Her ring|1963[|Schoeberl and Lindzen|1984 ). 



The EOMs for the two cumulants are integrated numerically using the same algorithms and 
methods as those employed for DNS, starting from the initial conditions ci(r) = Cjet(r) and 
c 2 (r, r') = c S 2 (r — r') — c/Atx with small positive c. The cumulants evolve toward the fixed 
point 

dci(r) = <9c 2 (r,r') = Q 
dt dt 

As a practical matter, we consider that the fixed point has been reached when the cumulants 
do not change significantly with further time evolution. It is essential for the second cumulant 
to have an initial non-zero value as otherwise it would be zero for all time, corresponding to 
axisymmetric flow, which is unstable with respect to non-axisymmetric perturbations. 

The programming task of solving the equations of motion for the cumulants is simplified by 
implementing the CE as a subclass of the DNS class, inheriting all of the lattice DNS methods 
without modification. The azimuthal symmetry of the statistics, Eqs. ( 14) and ( [T7j ), and the dis- 



crete symmetries, Eqs. ( [15] ) and ( ]18] ), are exploited to reduce the amount of memory required to 
store c 2 and p 2 . The symmetries also speed up the calculation and help thwart the development 
of numerical instabilities. The time step At is permitted to adapt, increasing as the fixed point is 
reached. Various consistency checks on the numerical solution are performed during the course 
of the time integration. For instance we verify that 

c 2 (r,r) = c 2 (0,0,AA = 0) > (25) 
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at all lattice points r. Furthermore from Eq. (fT2l) it must be the case that 



J Ci(r) dr = J C2(r, r') dr = 



(26) 



Likewise, as the second-order cumulant expansion conserves Kelvin's impulse, it follows from 



the impulse constraint (11) that 



Ci(r) sin(0) dr = q- ]C t{4>) sin(0) dr 



and 



y C2(r, r') sin(0) = . 

Finally the local mean kinetic energy is non-negative, 

(E(r)) = -p 2 (r,r) -j?i(r)ci(r) > 0, 



(27) 



(28) 



(29) 



because the statistics governed by Eqs. pTj ) and ( [22] ) are realizable: They can be obtained from a 
linear equation of motion for the vorticity that is driven by Gaussian stochastic forcing ( |Orszag 
1977 1 Salmon|1998 ). The total kinetic energy obtained by CE compares well to that determined 
by DNS. 



5. Comparison between DNS and CE 

The equal-time statistics accumulated in the DNS can be directly compared to the results of 
the CE because both calculations are based on the same jet model with the same finite-difference 
approximations on the same M x N lattice. Thus any differences between the DNS and CE 
statistics may be ascribed solely to the closure approximation. Results similar to those below 
are obtained on a coarser 200 x 100 lattice. 

Figure [5] shows the mean absolute vorticity calculated with the two approaches. Closest 
agreement between DNS and the CE is found at the shortest relaxation time of r = 1.5626 days. 
The CE is accurate for short relaxation times because fluctuations are suppressed by the strong 
coupling to the fixed jet. The second cumulant is reduced in size, and errors introduced by the 
closure approximation that neglects the third cumulant are small. For longer relaxation times, 
the CE systematically flattens out the mean absolute vorticity in the center of the jet too strongly. 
The largest absolute discrepancy in the mean vorticity appears at an intermediate relaxation time 
of r = 3.125 days. At longer relaxation times, the mean absolute vorticities in the DNS and CE 
become small in the central jet region; however, their fractional discrepancy increases, and the 
second cumulants show increasing quantitative and even qualitative discrepancies. 

Comparison of the second cumulants for r = 1.5625 days (Fig. [6]) reveals a qualitative 
discrepancy. The two-point correlations as calculated in the CE exhibit wavenumber-three 
periodicity, in disagreement with the wavenumber-four periodicity of the critical-layer wave 
dominating the fluctuating flow component in DNS (Fig [3]). In this regard, the CE mimics the 
wavenumber-three periodicity found in DNS at the longer relaxation times. In both DNS and 
CE, the correlations are strongest in absolute value when one of the two points of the second 
cumulant is located near the equator. Interestingly, the second cumulant from the DNS exhibits 
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a near-exact symmetry that is not a symmetry of the EOM, 



c 2 (-0,0',AA)^c 2 (0,0 / ,AA), (30) 



in addition to the symmetries of Eqs. ( 15 ) and ( 18 ). This approximate symmetry, which holds 



exactly for the second-order CE, may be attributed in the case of the DNS to the small size of 



the third cumulant. The fixed point of the second-order CE as described by Eqs. (21 ), (22 ), and 



( 24 ) possesses the artificial symmetry, for under the north-south reflection <p — > —<p the Jacobian 



operator ^ changes sign, as do both Ci((j>) and pi((f), and the fixed point equations remain 
unchanged provided that the second cumulant obeys Eq. ( |30| ). The artificial symmetry would, 
however, be broken in general by any coupling of the second cumulant to a third (non-zero) 
cumulant or, equivalently, by the inclusion of eddy-eddy interactions, which can redistribute 



eddy enstrophy spatially. Thus the artificial symmetry d30j) is an artifact of the closure (23 ). 



Other qualitative discrepancies appear at longer relaxation times (Fig. [7]). For r = 25 days, 



the second cumulant calculated by DNS no longer shows the artificial symmetry (30), whereas 
the symmetry continues to be present in the CE due to the closure approximation. In contrast to 
the r = 1.5625 case, the largest two-point correlations occur when one of the two points is away 
from the equator, reflecting the fact that correlations are washed out by the strong turbulence 
near the jet center. Finally, the second cumulant as calculated by CE shows a wavenumber-three 
periodicity, with excessively strong correlations at large separations, as a result of the neglect 
of eddy-eddy interactions, which strongly distort the eddy fields in the DNS. Nonetheless, even 
for relatively long relaxation times for which differences between the CE and the DNS at the 
center of the jet are apparent, the CE does capture the structure of the transition from the mixing 
region in the center of the jet to the non-mixing region away from the center, where the mean 
absolute vorticity in the DNS and the absolute vorticity of the underlying jet coincide. 

Figure [8] compares the eddy diffusivity k of absolute vorticity in the meridional direction 
as calculated by DNS and CE for the two cases of r = 3.125 and 25 days. The diffusivity is 
defined by second-order statistics: 

k(0 = -(v'(v)q'(v)) (DNS) 

l 

(CE). (31) 



d(a<p) 
1 dp 2 {4>, <f), A A) / dq 



cos0 <9(AA) 



The diffusivities calculated by the two methods are qualitatively similar: diffusion is largely 
confined to the mixing region of the jet; it is negative near critical layers within the mixing 
region (cf. Scho eberl andLindzen|1984| ). However, there are significant quantitative differences 



between DNS and CE, as might be expected given that second-order statistics will generally be 
poorly reproduced in a second-order closure. 



6. Discussion and conclusions 

The barotropic flows considered here attain statistically steady states after sufficient time 
has passed. They are out of equilibrium on large scales as the fixed zonal jet to which they 
relax is both a source and a sink of energy. Statistical approaches that have been developed to 
describe the equilibrium states of geophysical flows in the absence of large-scale forcing and 
dissipation therefore are not applicable here. For example, approaches based on maximizing 
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an entropy functional subject to constraints on energy, enstrophy, and possibly higher-order 
inviscid invariants (Miller] [T990| |Robert and Sommeria||1991[ |Salmon|[l998| |Turkington et al 



2001 1 Weichman~[2 006 ; Majda and Wang 2006) assume ergodic mixing and therefore would give 



statistical equilibrium states with mixing throughout the domain, rather than mixing confined 
to the region in the center of the jet. 

One-point closures likewise are generally not adequate for the inhomogeneous flows we 
considered. For example, as one ingredient of a diffusive one-point closure one might postulate 
a linear relationship between an eddy diffusivity k{4>) and rms fluc tuations in the streamfunction 



(ip 2 (r)} — (-^(r)) 2 , which has the same physical dimension (Holloway and Kristmannsson 



1984, Holloway|| 1986 [Stammer| 1998[ ). Figure [8] reveals that fluctuations in the streamfunction 



persist to much higher latitudes owing to Rossby wave dynamics outside the mixing region. A 
diffusive closure in which the diffusivity is a linear function of rms streamfunction fluctuations 
would therefore also lead to mean states with mixing in an overly large fraction of the domain. 

We have tested the simplest non-trivial two-point closure approximation based on a cumu- 
lant expansion of vorticity statistics, obtained by discarding the third and higher cumulants. 
This corresponds to neglecting eddy-eddy interactions while retaining eddy-mean flow interac- 
tions. The second-order CE is realizable, has no adjustable parameters (such as eddy damping 
terms), and is not restricted to homogeneous flows. For short relaxation times, the expansion 
reproduces the first moment fairly accurately. For longer relaxation times, it is quantitatively 
less accurate, but it still captures the transition from a mixing region at the center of the jet to a 
no-mixing region away from the center. 

The steady-state statistics from the CE can be found with much less computational effort 
than that required to calculate time-averaged statistics using DNS, as the partial differential 
equations governing the fixed point ([24]) are time-independent. This is especially true if a good 
initial guess is available for the cumulants c\ and c 2 because the fixed point can then be reached 
rapidly by iteration. Furthermore, as the statistics vary much more slowly in space than any 
given realization of the underlying dynamics (see Fig. [2]), it may be possible to employ coarser 
grids without sacrificing accuracy. Thus the CE realizes a program envisioned by Lorenz ( 1967 ) 
long ago by solving directly for the statistics, but it does so at the cost of a closure approximation 
that compromises the accuracy of the statistics, especially for flows with more strongly nonlin- 
ear eddy-eddy interactions. There is evidence, however, that nonlinear eddy-eddy interactions 
in Earth's atmosphere are weak ( |Schneider and Walk er|2006] ). |Q' Gorman and Schneider (2007) 
have shown that several features of atmospheric flows, such as scales of jets and the shape of the 
atmospheric turbulent kinetic energy spectrum, can be recovered in an idealized general circu- 
lation model in which eddy-eddy interactions are suppressed. In particular, baroclinic jets form 
spontaneously and are maintained by eddy-mean flow interactions even in the absence of non- 
linear eddy-eddy interactions. So a second-order CE may be worth exploring for more realistic 
models. Extensions of the CE to multi-layer models governed by either quasigeostrophic dy- 
namics or the primitive equations are straightforward. The storage requirement for the second 
cumulant grows as the square of the number of layers, but it remains feasible so long as the mod- 
els retain azimuthal symmetry. The incorporation of topography and other symmetry-breaking 
effects is more problematic because of the much larger storage requirement. 

Whether more sophisticated closures can be devised that are more accurate and yet only 
require comparable computational effort remains an open question. In the case of isotropic 
turbulence, renormalization- group inspired closures show some promise ( McComb~|2004[ ) , but 
these typically make extensive use of translation invariance in actual calculations. Investigation 
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of more sophisticated approximations for non-isotropic and inhomogeneous systems, such as 
the barotropic flows we considered, may be warranted in view of the partial success of the 
cumulant expansion reported here. 
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List of Figures 



Mean zonal component of the velocity (u((f))) on a unit sphere of radius a = 1. 
Relaxation times of r = [corresponding to the profile of the fixed jet Wj e t(0)], 

3.125 days, and 25 days are plotted [[6] 

Absolute vorticity q as calculated by DNS for a relaxation time of r = 25 days. 
The left and right hemispheres are shown in each panel; each is inclined by 
20° to make the poles visible. Deep red (blue) corresponds to q = ±1CT 4 s -1 . 
(a) Initial state with equatorial zonal jet. (b) Early development of instability, 
(c) Statistically steady state, (d) Mean absolute vorticity qi(r) = (q(r)) in 
statistically steady state, showing the effect of turbulence on the mean absolute 

vorticity profile and the recovery of azimuthal symmetry in the statistic [17] 

Snapshots of absolute vorticity in statistically steady states in a cylindrical pro- 
jection. The relaxation times are (a) r = 1.5625, (b) 3.125, (c) 6.25, (d) 12.5, 
(e) 25, and (f) 50 days. As in Fig. [I] deep red (blue) corresponds to q = ±10~ 4 



s- 1 . 



Different initial conditions yield the same low-order equal time statistics. The 
case of relaxation time r = 25 days is illustrated, (a) Lightly perturbed ini- 
tial absolute vorticity (from Fig. [T]). (b) Second cumulant obtained from the 
lightly perturbed initial condition with reference point (orange square) posi- 
tioned along the central meridian (A' = 0) and at latitude <p' = 18°. Colors 
indicate positive (deep red is 1CT 10 s~ 2 ) and negative (deep blue is — 1CT 10 s~ 2 ) 
correlations with respect to the reference point, (c) Highly perturbed initial 
condition, (d) Second cumulant obtained from the highly perturbed initial con- 
dition, (e) Comparison of the zonally averaged mean absolute vorticity in the 

central jet region [19] 

(a) Mean absolute vorticity q% as a function of latitude for different relaxation 
times. Zonally averaged results from DNS (solid lines) are compared to those 
from the CE (dashed lines). The black line (r = 0) is the absolute vorticity 
of the fixed jet gj et (0). (b) Magnified view of central jet region. Note the an- 
tisymmetry of the mean absolute vorticity (the first cumulant) under equatorial 

reflections l20l 

The second cumulant of the relative vorticity field, C2(0, 0', A — A'), for relax- 
ation time r = 1.5625 days, (a), (b) and (c): DNS. (d), (e), and (f): CE. The 
reference point (orange square) is positioned along the central meridian (A' = 0) 
and at latitudes of <p' = for (a) and (d), <p' = 18° for (b) and (e), and 4/ = 36° 
for (c) and (f). Colors indicate positive (deep red is 1CT 10 s~ 2 ) and negative 
(deep blue is — 1CT 10 s~ 2 ) correlations with respect to the reference point. . . . I2T1 
Same as Fig. [5] except for a relaxation time of r = 25 days. The reflection 
symmetry about the equator seen in the CE, an artifact of the closure truncation, 

is not present in the DNS [22] 

(a) Eddy diffusivity k for the case of relaxation time r = 3.125 days on the unit 
sphere (a = 1) as calculated by DNS and CE. Also plotted for comparison are 
the DNS calculation of the square root of the variance in the streamfunction, 
■>/ (ijj 2 (r)) — {ijj(r)) 2 , and the first moment of the absolute vorticity. (b) Same 
as (a) except for r = 25 days [23] 
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Fig. 1. Mean zonal component of the velocity (u((f>)) on a unit sphere of radius a — 1. Relax- 
ation times of r = [corresponding to the profile of the fixed jet M jct (0)], 3.125 days, and 25 
days are plotted. 
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Fig. 2. Absolute vorticity q as calculated by DNS for a relaxation time of r = 25 days. 
The left and right hemispheres are shown in each panel; each is inclined by 20° to make the 
poles visible. Deep red (blue) corresponds to q = ±1CT 4 s^ 1 . (a) Initial state with equatorial 
zonal jet. (b) Early development of instability, (c) Statistically steady state, (d) Mean absolute 
vorticity gi(r) = (q(r)) in statistically steady state, showing the effect of turbulence on the 
mean absolute vorticity profile and the recovery of azimuthal symmetry in the statistic. 



17 




ISO 120W 60W 60E 120E 180 ISO 120W 60W 60E 120E 180 




180 120W 60W 60E 120E 180 1 80 120W AOW W1E 120E 180 




180 120W GOW 60E 120E 180 180 120W 60W 60E 120E 180 



Fig. 3. Snapshots of absolute vorticity in statistically steady states in a cylindrical projection. 
The relaxation times are (a) r = 1.5625, (b) 3.125, (c) 6.25, (d) 12.5, (e) 25, and (f) 50 days. 
As in Fig. [TJ deep red (blue) corresponds to q = ±10~ 4 s _1 . 



18 




180 120W 60W 60E 120E 180 180 120W 60W 60E 120E 180 







-6 



; (e) 






































































































Lightly Perturbed DNS 










-•- Highly Perturbed DNS 


























25 -20 -1 


5 -1 





5 E 


1 


1 


5 20 2 



Latitude (degrees) 

Fig. 4. Different initial conditions yield the same low-order equal time statistics. The case 
of relaxation time r = 25 days is illustrated, (a) Lightly perturbed initial absolute vorticity 
(from Fig. [T]). (b) Second cumulant obtained from the lightly perturbed initial condition with 
reference point (orange square) positioned along the central meridian (A' = 0) and at latitude 
</>' = 18°. Colors indicate positive (deep red is 1CT 10 s~ 2 ) and negative (deep blue is — 10~ 10 
s~ 2 ) correlations with respect to the reference point, (c) Highly perturbed initial condition, (d) 
Second cumulant obtained from the highly perturbed initial condition, (e) Comparison of the 
zonally averaged mean absolute vorticity in the central jet region. 
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Fig. 5. (a) Mean absolute vorticity q\ as a function of latitude for different relaxation times. 
Zonally averaged results from DNS (solid lines) are compared to those from the CE (dashed 
lines). The black line (r = 0) is the absolute vorticity of the fixed jet 9j et (</>). (b) Magnified view 
of central jet region. Note the antisymmetry of the mean absolute vorticity (the first cumulant) 
under equatorial reflections. 



20 




60N- 
30N- 
EQ- 
30S 



(C) 




ii I20W 


60W 


60E I20E ) 


(e) 








I20W 60W 60E 120E 1 


(f) 




(1 120W 


60W 





Fig. 6. The second cumulant of the relative vorticity field, c 2 (0, 0', A — A'), for relaxation time 
r = 1.5625 days, (a), (b) and (c): DNS. (d), (e), and (f): CE. The reference point (orange 
square) is positioned along the central meridian (A' = 0) and at latitudes of <// = for (a) and 
(d), 0' = 18° for (b) and (e), and 4>' = 36° for (c) and (f). Colors indicate positive (deep red 
is 10"" 
point. 
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Fig. 8. (a) Eddy diffusivity k for the case of relaxation time r = 3.125 days on the unit sphere 
(a = 1) as calculated by DNS and CE. Also plotted for comparison are the DNS calculation of 
the square root of the variance in the streamfunction, a/ (^ 2 (r)) — (^(r)} 2 , and the first moment 
of the absolute vorticity. (b) Same as (a) except for r = 25 days. 
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